<< Interplanetary cruise Cookbook Coordinates and frames >>

CelestLab >> - Introduction - > Cookbook > Jacobian and covariance matrices

Jacobian and covariance matrices

Jacobian and covariance matrices

Introduction to the examples

The following examples show the use of jacobian and covariance matrices in concrete applications.

Quick recap:

The following examples show:

Note: First execute the code below to initialize some variables used in the examples.

// -----------------------------------------------
// Data
// -----------------------------------------------

// Initial date
cjd0 = CL_dat_cal2cjd(2015,1,1);

// Associated orbital elements (Keplerian): semi-major axis, eccentricity, 
// inclination, argument of perigee, RAAN, mean anomaly (units = m, rad)
sma = 7.e6; 
ecc = 0.001;
inc = CL_deg2rad(45);
pom = 0;
gom = 0;
anm = 0;
kep0 = [sma; ecc; inc; pom; gom; anm]; 

// Propagation dates (2 days)
cjd = cjd0 + (0 : 1/1440 : 1); 

// Ground station (geodetic coordinates)
sta = [0; CL_deg2rad(45); 0]; 

// Station cartesian coordinates
psta = CL_co_ell2car(sta); 

// Station local vertical
zsta = CL_co_sph2car([sta(1:2); 1]);

// Minimum elevation
elevmin = CL_deg2rad(5);

Examples

1) Basic use of Jacobian matrices

The objective here is to compute the Jacobian of the transformation from orbital elements to (inertial) position and velocity in the local frame.

// -----------------------------------------------
// Keplerian elements to local frame
// -----------------------------------------------
kep = kep0; 

// kep -> pv
[pos, vel, dpv_dkep] = CL_oe_kep2car(kep); 

// pv -> pv in local frame
M = CL_fr_tnwMat(pos, vel); 
J =  [M, zeros(M); zeros(M), M]; 

// Jacobian of combined transformation
dpvloc_dkep = J * dpv_dkep 

// Note that the result could also be obtained by: 
dpvloc_dkep2 = CL_op_orbGapLofMat("kep", kep, "tnw", meth="c", form="un")

2) Jacobian matrices in conjunction with orbit propagation

Same as example 1 except that the Keplerian elements are obtained by propagation.

// -----------------------------------------------
// Keplerian elements to local frame (propagation)
// -----------------------------------------------

// kep0 -> kep
[kep, dkep_dkep0] = CL_ex_secularJ2(cjd0, kep0, cjd); 

// kep -> pv 
[pv, dpv_dkep] = CL_oe_convert("kep", "pv", kep); 

// Jacobian of combined transformation
dpv_dkep0 = dpv_dkep * dkep_dkep0; 

// Extract d(pos) / d(sma, inc)
dpos_dai = dpv_dkep0(1:3,[1,3],:); 

// Inertial frame -> local frame
M = CL_fr_tnwMat(pv(1:3,:), pv(4:6,:)); 

// Jacobian of transformation: (sma, inc) -> position in local frame
dpvloc_dai = M * dpos_dai; 

// Plot sensitivity of position in local frame to semi major axis (unit: m/m)
scf(); 
plot(cjd-cjd0, matrix(dpvloc_dai(:,1,:),3,-1)); 
CL_g_stdaxes(); 
CL_g_legend(gca(), ["t", "n", "w"]); 

// Plot sensitivity of position in local frame to inclination (unit: m/m)
scf(); 
plot(cjd-cjd0, matrix(dpvloc_dai(:,2,:),3,-1) / sma); 
CL_g_stdaxes(); 
CL_g_legend(gca(), ["t", "n", "w"]);

3) Basic use of covariance matrices

Starting from example 1, the covariance matrix for the position in the local frame is expressed as a function of the covariance matrix for the orbital elements.

// -----------------------------------------------
// Keplerian elements to local frame (covariance matrix)
// -----------------------------------------------

// Same as example 1 - see above for explanations
kep = kep0; 
[pos, vel, dpv_dkep] = CL_oe_kep2car(kep); 
M = CL_fr_tnwMat(pos, vel); 
J =  [M, zeros(M); zeros(M), M]; 
dpvloc_dkep = J * dpv_dkep; // size = 6x6

// Covariance matrix for the orbital elements 
// 1-sigma errors (arbitrary values)
sig_sma = 1000; // m
sig_ecc = 1000 / sma; 
sig_inc = CL_deg2rad(0.01); // rad 
sig_pom = CL_deg2rad(0.01); 
sig_gom = CL_deg2rad(0.01); 
sig_anm = CL_deg2rad(0.1); // <=> ~12 km (tangential axis)

// Diagonal covariance matrix
cov_kep = diag([sig_sma, sig_ecc, sig_inc, sig_pom, sig_gom, sig_anm].^2) 

// Covariance of position/velocity errors in local frame
cov_pvloc = dpvloc_dkep * cov_kep * dpvloc_dkep' 

// Standard deviations of position errors for each axis (t, n, w) 
sig_errpos = diag(sqrt(cov_pvloc(1:3, 1:3))) // meters

4) More advanced use of covariance matrices

The objective is to evaluate the acquisition angular tolerance by a ground station after launch.

Only injection errors on semi-major axis and inclination ar supposed to have an sensitive impact.

// -----------------------------------------------
// Acquisition tolerance evaluation
// -----------------------------------------------
kep = kep0; 

// First propagate the orbit and compute visibility intervals 
[kep] = CL_ex_secularJ2(cjd0, kep0, cjd); 
[pos, vel] = CL_oe_kep2car(kep); 
posecf = CL_fr_convert("ECI", "ECF", cjd, pos, vel); 
intvisi = CL_ev_stationVisibility(cjd, posecf, sta, elevmin); 

// Consider beginning of intervals (potential acquisition dates) 
cjdv = intvisi(1,:); 
N = size(cjdv, "*"); 

// Compute impact of initial semi-major axis and inclination 
// on position perpendicular to line of sight

// Propagation
[kep, dkepdkep0] = CL_ex_secularJ2(cjd0, kep0, cjdv); 

// Cartesian coordinate (ECI then ECF)
[pos, vel, dpvdkep] = CL_oe_kep2car(kep); 
[posecf, velecf, dpvecfdpv] = CL_fr_convert("ECI", "ECF", cjdv, pos, vel); 

// Station local frame (Z axis = line of sight, X axis = local vertical) 
M = CL_rot_defFrameVec(posecf-repmat(psta,1,N), zsta, 3, 1); 

// Jacobian of transformation: (sma0, inc0) -> perpendicular position error
// => J: 2x2xN matrix
J = M(1:2,:,:) * dpvecfdpv(1:3,:,:) * dpvdkep * dkepdkep0(:,[1,3],:); 

// Assumed errors on initial sma and inc
covai = diag(([1.e3, CL_deg2rad(0.1)]).^2); 

// Covariance matrix on (x,y) perpendicular to line of sight
// Note: covai is duplicated N times (=> 2x2xN)
covacq = J * repmat(covai, [1,1,N]) * J'; 

// Show angular errors in degrees 
for k = 1 : N
  covang = ((covacq(1:2,1:2,k) ./ (CL_norm(posecf(:,k) - psta))^2)); 
  ang = atan(sqrt(max(real(spec(covang))))); 
  mprintf("%d: %f %f\n", k, cjdv(k) - cjd0, CL_rad2deg(ang));  
end


Report an issue
<< Interplanetary cruise Cookbook Coordinates and frames >>